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ABSTRACT 


The study describes two problems related to radar rain 
fall estimation. The first part is a description of a pre- 
liminary data analysis for the purpose of statistical estima- 
tion of rainfall from multiple (radar and raingage) sensors. 
Raingage, radar and joint radar-raingage estimation is de- 
scribed and some results are give. Statistical parameters of 
rainfall spatial dependence are calculated and discussed in 
the context of optimal estimation. Quality control of radar 
data is described also. The second part describes radar 
scattering by ellipsoidal raindrops. Analytical solution is 
derived for the Rayleigh scattering regime. Single and vol- 
ume scattering is presented. Comparison calculations with 
the known results for spheres and oblate spheroids are shown. 
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PART I 


ESTIMATION OF RAINFALL FROM RADAR AND RAINGAGES 
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1.1. INTRODUCTION 


Estimation of global rainfall is an important component 
of global climate studies. This has been well established 
and documented (Simpson [1989], Thiele [1987], Wilkerson 
[1988], Arkin and Ardunay [1990]) . Due to the fact that over 
70% of the Earth's surface is covered by oceans it is neces- 
sary to use satellite technology for global rainfall estima— 
tion. Satellite methods of rainfall estimation rely on indi- 
rect ways of inferring rainfall over an area based on mea- 
surements of radiation emmitted from several different fre- 
quency bands . For review of satellite methods of rainfall 
estimation refer to Barrett and Martin [1981] , Adler and 
Negri [1988], or Arkin and Ardunay [1990]. Satellite estima- 
tion methods require validation because they are based on in- 
direct inference. Validation simply means comparing the es 
timation with other independently obtained estimates of the 
same rainfall. The validation will not be meaningful if the 
reference methods are of poor quality or if their accuracy 
cannot be established. The reference methods are typically 
based on raingage networks and more recently on estimates 
from weather radars. It is the purpose of this report to 
discuss these methods with particular emphasis on two as- 
pects: 
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1 . 


estimation of mean areal rainfall using a combina- 
tion of radar and raingage observations; and 

2. estimation of uncertainty associated with areal es- 
timates of rainfall 

A general discussion of the problem will be illustrated 
with analysis of radar rainfall data from Darwin, Australia. 
The data set used in the study is fully described in 
Krajewski and Rexroth [1990] . 

1.2. RAINGAGE RAINFALL ESTIMATION 

Raingages are the most traditional means of measuring 
rainfall. The observed rainfall represents a point value and 
in order to make an assessment of areal rainfall one needs to 
resort to interpolation methods. Mathematically the problem 
can be defined in the following way. Suppose that the true 
but unknown mean areal rainfall is 

R a = 7 1 p (u, t ) du (1) 

where A is the area of interest, u=(x,y) is a two dimensional 
space location, and P is the true rainfall process. The in- 
dex t signifies the temporal aspect of the problem. In the 


3 


following derivations this index will be dropped to simplify 
the notation. It will be assumed that both the estimated 
areal rainfall and its observations are given at the same 
time scale. The temporal resolution of the data used in this 
work is either one hour or one day. Therefore, equation (1) 
denotes, say, daily mean area rainfall. The task at hand is 
to estimate R A given n point raingage observations Gi(u,t), 

i=l,2,...,n. A linear estimator of the following form can be 
proposed 

ft A = I XiGi(u,t) (2) 

i-1 

It is required that the estimator ft A be unbiased and give 
minimum variance. The unbiasedness condition gives 

E [ £ XiGi(u,t) ] - E [J f P(u,t)du] (3) 

i-l A 

By taking advantage of the linearity of both sides of this 
condition it can be simplified to 

£ XjE [Gi(u,t) ] J E [p (u, t) du] (4) 

i-l A 


Assuming that the mean of the rainfall process P is constant 
in space and that the observations Gi(u,t) are unbiased, 

equation (4) can be reduced to 
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( 5 ) 


n 


I Xi - 

i-1 


1 


Minimization of the variance of the estimator is accom- 
plished in the following way. First, the variance of the es- 
timator is written as 

Or - E[(R A -ft A ) 2 ] = E[R A 2] - 2E[R a & a ] + E[ft A 2 ] 

= E[A- J / p (Ui) )P (u 2 )du 1 du 2 ] 

A A 
1 ^ 

- 2E[7 J X XiGi(u)P (u)du] 

A a i-i 

+ E[ £ X XiXjGi (u) G j (u) ] ( 6 ) 

i-1 i-1 


Making use of the fact that for the second-order stationary 
processes 


cov(v) - E[P(u 1 )P(u 2 ) ] - m 2 = cov(u 1 -u 2 ) 


where m is the mean of the process, equation (6) can be writ 
ten as 


2 1 

Or = ~2 J J cov (u^— u 2 ) du x du 2 

A A A A 


- \ J X XiCOvtu-Ui) du 
A A i-1 
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( 8 ) 


n n 

+ £ £ X i X j cov(u i -Uj) 

i~l i*l 

The next step is finding the weights that will minimize 
the expression (8) subject to the constraint (5) . This could 
be accomplished using the method of Lagrange multipliers. 
According to the method an unconst raind optimization problem 
is solved which minimizes the following function 

F = Or + 2|l( £ A* - 1) (9 ) 

i-l 


This leads to the following set of linear equations 


n 


= J cov(u-Ui)du + 2 £ XjCOv(Ui-Uj) + 2\l = 0 
dXi A a j-1 

for i=l, 2, ... ,n 


( 10 ) 


-£•- £ n, - 1 

3Xi "i 


(ID 


The solution of this system yields the optimal set of 

★ , 
weights denoted with X i which when substituted into the vari- 
ance equation (8) give 

1 n * 

= ~^2 \ i cov (U!“U 2 ) dUidU 2 - ^ J £ A^covfu-Uj.) du - \l (12) 
A A A A A i-l 
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In order to use the above method one needs to estimate 
the spatial covariance function of the rainfall process. Of 
course the true covariance, which appears in all the double 
integral terms, cannot be inferred from the noisy data and 
w ill be substituted by a model obtained by fitting a theoret- 
ical covariance function to the data. 

There are many functions which could be used as models 
for covariance. For a discussion of these models and the 
conditions they have to meet see Journel and Huijbregts 
[1978] . Before a model is decided on, the empirical (or raw) 
covariance values need to be inspected. In Figures 1 and 2 
the raw correlations (normalized covariances) are shown for 
all the pairs of raingages for daily and hourly data, respec- 
tively. The network of raingages under Darwin radar is com- 
posed of 22 raingages (21 were included in the analysis) . 
Their locations and the raingage data are discussed in 
Krajewski and Rexroth [1990] . It is clear from the Figures 
that there is not much correlation between the rainfall val- 
ues recorded at the network gages. There are two possible 
explanations. The first possibility is the statistical sam- 
pling variability. The data used in the analysis cover a 
rather 
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1. Experimental correlation obtained from daily 

raingage data for the period December 21, 1987 

- January 20, 1988. 
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Experimental correlation obtained from hourly 
raingage data for the period December 21, 1987 
- January 20, 1988 
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CORRELATION 



DISTANCE (km) 

Figure 3- Experimental correlation obtained from daily 

raingage data for the period October 31, 1987 - 
May 15, 1988. 
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short period of time (1 month, December 21, 1987 through 

January 20, 1988) . It is then conceivable that the lack of 

significant correlation can be attributed to the small sample 
size used. To verify this hypothesis correlation was calcu- 
lated for the 6 months of data available covering the period 
from October 31, 1987 to May 15, 1988. The results are pre- 

sented in Figure 3. The second possibility is that the ex- 
isting network is to sparse to capture the fluctuation scale 
of the events. The shortest distance between any two sta- 
tions within the network is 14 km and there are only a few 
pairs of stations separated by a distance less that 30 km. 
Correlation distance defined as distance at which the expo- 
nential correlation drops to l/e=0.37 depends on the temporal 
scale of interest. For hourly GATE data the correlation dis- 
tance is about 20 km (Bell [1987]). For daily GATE data it 
increases to about 40 km. These numbers represent an overes- 
timation of the point value correlation since they were cal- 
culated from radar data averaged in space to 4 by 4 km val- 
ues. Correlations were also calculated from the Darwin 
radar-rainfall data. They are shown in Figures 4 and 5. As 
a conclusion from the above plots and the discussion it is 
fair to suspect that it is the sparseness of the network that 
causes the lack of significant correlation. It can be shown 


that, if there is no significant correlation the weights are 
all equal and the estimator becomes a simple average. 
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Figure 4 . 


Experimental correlation obtained from hourly 
radar data for the period December 21, 1987 - 
January 20, 1988. 



Y-Correlation 



Figure 5 . 


Experimental correlation obtained from daily 
radar data for the period December 21, 1987 

January 20, 1988. 
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X-Correlation 


Based on the raingage data the average monthly rainfall 
rate for the Darwin radar umbrella, calculated from the rain- 
gage data is 0.45 mm/hr. The values of daily rainfall are 
given in Table 1. The variances of these rainfall estimates 
can be calculated simply as 

Var {R a -P } = 1/n [<T 2 -E {Cov (X, Y) }] (13) 

where X,Y are independent uniformly distributed points 
(gages) in the domain of interest. If this domain is large 
compared to the effective range of the covariance function 
the second term in (13) becomes negligible. In our case the 
domain of interest has radius of 170 km while the effective 
range of the correlation is only about 20 km. Assuming tem- 
poral independence of the daily estimates the standard devia- 
tion of the monthly estimate of the average rainfall rate is 
about 0.02 mm/hr. This corresponds to about 13 mm for 
monthly rainfall accumulation. 
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Table 1. Daily Areal Averages of Radar and Gauges 


Date 


R-Rate 

(mm/hr) 


R-Acum 

(mm) 


G-Rate 

(mm/hr) 


G-Acum 

(mm) 



12/21/87 

0.8937 

12/22/87 

1.0331 

12/23/87 

1.1231 

12/24/87 

0.5020 

12/29/87 

0.2238 

12/30/87 

0.4721 

12/31/87 

0.9245 

01/01/88 

0.0062 

01/02/88 

0.2696 

01/03/88 

0.2572 

01/04/88 

0.3915 

01/05/88 

0.2628 

01/06/88 

0.1724 

01/07/88 

0.1262 

01/08/88 

0.1488 

01/09/88 

0.3547 

01/10/88 

0.5545 

01/11/88 

0.5127 

01/12/88 

0.1859 

01/13/88 

0.1783 

01/14/88 

0.1640 

01/15/88 

0.0867 

01/16/88 

0.3599 

01/17/88 

0.0752 

01/18/88 

0.0484 

01/19/88 

0.2773 

01/20/88 

0.2725 


10.7 1.1823 28.4 

7.2 1.3912 33.4 

20.2 1.5884 38.1 

2.0 0.1164 2.8 

2.5 0.0794 1-9 

11.3 1.0933 26.2 

22.2 1.4980 36.0 

0.1 0.0000 0.0 

6.5 0.4370 10.5 

4.4 0.1951 4.7 

9.4 0.3968 9.5 

6.3 0.2799 6.7 

3.3 0.3857 9.3 

3.0 0.0590 1.4 

3.6 0.2360 5.7 

7.1 0.7094 17.0 

11.6 0.4539 10.9 

12.3 0.4026 9.7 

4.5 0.0919 2.2 

4.3 0.3670 8.8 

3.6 0.2223 5.3 

2.1 0.0097 0.2 

8.6 0.3957 9.5 

1.8 0.0164 0.4 

1.2 0.0000 0.0 

5.8 0.0919 2.2 

2.5 0.1092 2.6 


The monthly areal averages for radar and gages are : as *ol- 
lows* R-Rate-0.36 mm/hr, R-Accum=147 mm, G Rate 0.44 
G-Accum- 291 mm. The gage rates were calculated with the as 
sumption of 24 hr measurement by all ^ges during the entire 
period. This was not always the case. Missing hours were 
accounted for in the monthly and daily radar averages. The 
Channel gage was not included in the gage averages. 
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1.3. RADAR-RAINFALL ESTIMATION 


Estimation of Darwin rainfall based on the radar data 
was described in Krajewski and Rexroth (1990) . In their at 
las of radar— rainfall no attempt was made to use the best (in 
some sense) Z-R relationship or special preprocessing of the 
radar data. The 1.5 km CAP I data prepared at the NASA 
Goddard Space Flight Center were used. The Z-R relationship 
used was the following: a=230 and b=1.25. 

It was found that anomalous propagation presented a se- 
vere problem and therefore, a manual computer graphics-aided 
procedure was utilized to eliminated it as much as possible. 
Also, a quality control procedure aimed at detecting isolated 
outliers (see Krajewski, 1987) was applied to both hourly and 
daily rainfall fields converged to a 4 km by 4 km resolution 
grid. Outliers were defined as those data points which were 
statistically inconsistent with their immediate neighbor- 
hoods. The outliers detected in the daily fields are listed 
in Table 2. The critical parameter a which appears in the 

table controls the sensitivity of the method. Low values of 
the parameter indicate high sensitivity and, as a result, 
many points are questioned as being outliers. High values of 
a correspond to lower sensitivity of the method, and as a re- 
sult some of the "not-so-obvious" outliers may slip through 
the quality control. Based on the simulation study performed 
by Krajewski (1987) the optimal choice of the critical param- 
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eter for daily rainfall fields is a-3.0. The total number of 
outliers detected with this value of a in the studied period 

was 8 . 

Once outliers are detected the problem becomes how to 
accommodate them. By accommodation is meant replacement of 
the erroneous values with their surrogates. One "safe" 
choice of such replacement is to use the local mean of the 

surrounding data. 

After the quality control steps were implemented monthly 
rainfall based on radar observations was calculated from the 
daily fields. Comparison of daily and monthly values of 
radar-rainfall is given in Table 2. 
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Critical Parameter a 





Table 2. List of outliers detected in the daily data. 
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1.4. RADAR-RAINGAGE COMBINATION 


The time-area average of interest can be expressed math- 
ematically according to eguation (1) • An estimate of Rta 
noted Rta can be obtained from a number of ground-based sen- 
sors. The most appropriate network for the purpose of satel- 
lite rainfall validation exercise seems to be a combined net- 
work of a radar and a number of raingages. We could write 

R ta - g < 14 > 

where £r= (Zri, Zr2? • • • / Zrnr) a vector °f radar observa- 

tions within the area of interest A, .Zg* (Zgi, Zg2/ • • • > Zgng) 
a vector of NG corresponding raingage observations, and g is 
a function of NG and NR arguments . It is often taken to be 
linear, however it can also be nonlinear . The error term £ 

resulting from the approximation (14) is a random component 

2 

which could be characterized by its mean Jig and variance d £ . 
Below we outline a methodology to compute Rta an d the associ- 
ated Og from radar and raingage observations. This methodol- 
ogy is based on the main assumption that rainfall is a real- 
ization of a spatial stochastic process (random field) . 

Weather radars could be classified as indirect rainfall 
measuring devices in the sense that they measure parameters 
related to rainfall and not the rainfall itself. Rainfall 
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rate, and more precisely rainfall volume is, therefore, esti 
mated from those measurable parameters . As such, radar - rain- 
fall estimates are contaminated by errors which have two ba- 
sic components: 1) measurement error — error in measuring 

the parameters related to rainfall; and 2) inference error 
error due to an imperfect model (i.e., relationship) assumed 
between the measured parameters and rainfall. For example, 
one of the major errors introduced in inference models is ice 
contamination. The chances for ice contamination of the 
radar signal increase with increasing range. Contamination 
by ice increases reflectivity for the same rainfall, some- 
times by several dB . Another very important measurement er- 
ror that is normally neglected, but of extreme importance is 
radar calibration error. This is a systematic error unlike 
many other errors that are random in nature. This calibra- 
tion error can easily be of the order of couple dB, which re- 
sults in high errors in rainfall estimates. This radar cali- 
bration error can also drift with time, potentially changing 
with season. Figure 6 demonstrates the behavior of the bias 
which results from many different sources, not necessarily 
just the lack of calibration. Appendix C includes a simula- 
tion study of a bias model . (Krajewski and Smith, 1991). 
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G-Accum/R-Accum 


Figure 6 . 



Time series of daily bias of the radar-rainfall 


ThP bias is defined as the ratio of gage to 


rouA 


mean areal values . 
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Thus, radars which provide nearly continuous spatial 
coverage of large areas, are characterized by often signifi- 
cant point errors . On the other hand, a raingage network of- 
fers good point accuracy but is characterized by relatively 
high sampling error related to network density. The premise 
of a multisensor rainfall estimation is a combination of high 
spatial resolution of radar and good point accuracy of rain- 
gages . 

II. 4.1 Methodology outline 

The characteristics of radar and raingage sensors have 
been discussed in literature many times (see Wilson and 
Brandes, [1979]; Doviak and Zrnic, [1984]; Hudlow et al., 
[1984]; Zawadzki, [1982]; Austin, [1987]; Krajewski, [1987]) 
justifying the concept of a combined network. Several ap- 
proaches have been proposed to combine the two types of mea- 
surements (Brandes, [1975]; Crawford, [1979]; Eddy, [1979]; 
Krajewski, [1987]; Seo et al., [1990ab]; Azimi-Zonooz et al., 
[1988]; Smith et al., [1990]; and Seo and Smith [1990]). The 
recent works offer a comprehensive approach to the sensor 
merging problem. This approach is based on the stochastic 
interpolation technique called cokriging (Journel and 
Hui jbregts, [1978]). Its application requires estimation of 
spatial covariance functions of data from the involved sen- 
sors and the cross-covariance function. The method is capa- 
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ble of accounting for different sampling geometries of the 
measurements involved and the measurement error. It is as- 
sumed that 


NR NG 

g (2r, ZkG) = £ pRi z Ri + 21 PGj z Gj 

i=l 


(15) 


where the coefficients Pri and Poj can be found by minimizing 
the estimation error variance 


A A 

Var [Rta ~ r ta1 = E ( [ r ta - r ta1 


2 \ = 


(16) 


subject to the condition that we consider only unbiased esti 
mators, i.e. 


E [Rta! = r ta (17) 

The solution of this minimization problem is a function of: 
1) observations Z Ri and Z Gj ; 2) the spatial covariance func- 
tions of radar observations and raingage observations; 3) the 
spatial cross covariance function between the radar and rain- 
gage observations; and 4) the spatial covariance functions 
between the observations and the true rainfall. This last 
term is of course unknown and in general cannot be estimated 
from the data. Krajewski [1987] proposed a simple parameter- 
ization of this covariance term and performed a limited sen- 
sitivity analysis of such an approach. The results showed 
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quite flat behavior of the criterion function near the opti 
mal location of the parameters. However, the disadvantage of 
this approach is the lack of meaningful interpretation of the 
parameters, which in turn may cause difficulties in choosing 
the parameters for real- data situations . The covariance be- 
tween the observations and the true rainfall be expressed as 
a function of measurement error parameters of the two sen- 
sors. Radar observations can be expressed as 

T 

ZRi=ij 757 J R(t,s)dtds+coi i=l,...,NR (18) 

T 0 1 b 'S 

where S is the area of a basic radar observation (typically 
about 4 km x 4 km grid) , and COi is an error term. This error 

term can be characterized by its second order moments: the 

mean and the covariance. 

The raingage data represent point observations 
1 T 

Z G j=^J R (t , S) dt + vj j=l,...,NG (19) 

0 

2 

where vj is a zero mean Gaussian variable with variance <J V 

and can be assumed independent of R(t,s). Therefore, the co- 

variance between the observations and the true rainfall can 

2 

be expressed as a function of |lci)r covqj and O v . These parame- 
ters, in addition to the raingage network density and the 
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methodology used for combining the two data sets, control the 
performance of the system. 

The precise knowledge of these parameters results in an 
optimal scheme of merging radar and raingage data. However, 
in general it is very difficult to obtain the true values of 
these parameters. Therefore, it is important to investigate 
the sensitivity of the merging scheme with respect to the un- 
certainty about the measurement error parameters. Such a 
sensitivity study can be achieved by a simulation experiment. 

It is clear from the previous discussion that the error 
term of the estimate (14) could be directly used in the vali- 
dation of the satellite-based rainfall estimation methods. 

The problems of estimating the covariance function from 
a limited sample is demonstrated in Figures 7 and 8 and Table 
3. They show the range dependence of the cross-correalation 
function between the radar and raingage data. A strong ef- 
fect is evident. The cause was probably the fact that the 
Darwin radar uses a 5 cm, and therefore attenuating, wave- 
length . 
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Radar/Gage Correlation 



2 ' 



Table 3 . Daily Gage Correlations 


Gage # 

Gage Name 

Sanaa. 

Correlation 

Data Points 

1 

Annabur 

96.25 

0.724 

23 

2 

Bachelor 

67.27 

0.425 

22 

3 

Bathurst Island 

83.22 

0.521 

24 

4 

Bellville Park 

34.13 

0.733 

23 

5 

Charles Point 

31.58 

0.858 

24 

6 

Dum In Mirrie 

62.30 

0.789 

21 

7 

Garden Point 

125.72 

0.387 

20 

8 

Goodall Mine 

98.27 

0.747 

19 

9 

Gunn Point 

27.29 

0.937 

24 

10 

Humpty Doo 

44.78 

0.943 

25 

11 

Label le 

29.61 

0.842 

24 

12 

Koolpinyah 

89.87 

0.442 

23 

13 

Litchfield 

116.30 

0.292 

23 

14 

Mandorah Jetty 

17 . 12 

0.807 

23 

15 

McMinns Lagoon 

21.47 

0.867 

25 

16 

Mt . Bundey 

56.82 

0.292 

26 

17 

Old Point Stuart 

97.51 

0.675 

23 

18 

Pickertaramoor 

76.16 

0.929 

22 

19 

Point Stuart 

92.07 

0.723 

24 

20 

Snake Bay 

114.76 

0.282 

22 

21 

Woolner 

63.79 

0.781 

26 
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It should be evident from the shape of the correlation func- 
tion that as the separation lag increases the correlation at 
greater distances from the radar will decrease even faster. 
Therefore, the radius of the statistical influence of the 
raingage points on the radar data in the vicinity of the 
gages is relatively small. It should be expected that the 
radar-based estimates of rainfall will not be strongly af- 
fected by the gage estimates . 


1.5. CONCLUSIONS 


The performed analysis of the Darwin rainfall data for 
one monthly period demonstrated that a statistical approach 
to rainfall estimation using both radar and raingage data is 
faced with major problems. The two biggest problems seem to 
be : 

1. Quality control and preprocessing of radar data. 
Efficient anomalous propagation detection and re- 
moval is fundamental to the success of any subse- 
quent estimation approach. 

2. Sparseness of the raingage network and its configu- 
ration. Since the statistical estimation approach 
relies on the inference of statistical moments from 
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the observations, the sampling density is a criti- 
cal issue. The network of 22 approximately evenly 
distributed gages cannot capture the spatial vari- 
ability of the convective and monsoonal rainfall. 
However, the observed difficulty and apparent lack 
of strong statistical relation between the radar 
and the raingage data should be investigated fur- 
ther with a longer data set . The way the network 
configuration comes into play is through its abil- 
ity to observe the spatial correlation structure of 
the investigated process. In the case of the 
Darwin network it seems that the rainfall process 
has correlation distance on the order of 20 km or 
less even on the daily scale. The network which 
has average intergage separation distance greater 
then that cannot resolve such a scale. 

Investigation of the statistical approach as a nonsta- 
tionary and space— time process should and will continue . 
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PART II 


RADAR BACKSCATTERING BY ELLIPSOIDAL RAIN DROPS 
ANALYTICAL SOLUTIONS FOR RALYLEIGH REGION 


41 


II.l. INTRODUCTION 


The calculations of scattering of electromagnetic waves 
a.re required in many fields such as remote sensing or radia 
tive heat transfer. In particular, investigation of the 
scattering of radar waves by raindrops and the scattering of 
light by small chemical and biological particles is of inter- 
est in studies related to global climate modeling . 

The solution of the electromagnetic scattering problem 
by spherical objects is well known as the Mie theory (see 
[11,26,39]) and has been used to a great advantage in many 
physical applications. However, the problem of scattering by 
nonspherical bodies often arises. Many techniques have been 
developed but each has a limited range of applicability that 
is determined by the size of the scattering object relative 
to the wavelength of the incident field. 

The scattering by objects that are very small compared 
to the wavelength can be analyzed by Rayleigh approximation. 
Since Rayleigh's classical paper [35], there have been many 
studies performed in this area. Let us only mention 
Kleinman's work [29,30] for general considerations of the 
low-frequency electromagnetic scattering problems, Jones [25] 
for the theoretical aspects of the scattering problem, 
Kleinman and Senior [30] for their investigations on scatter- 
ing cross-sections, Asvestas and Kleinman [5,6] for the solu- 
tion of the low-frequency scattering problem by spheroids and 
disks, Angell and Kleinmann [1] for polarizability tensors, 
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Stevenson [37,38] for the solution in the case of the ellip- 
soid, Siegel [36] for work on bodies of revolution, and 
Darling and Senior [16] for a general consideration of low- 
f reguency scattering by separable and nonseparable bodies . 

Objects whose size is of the order of the wavelength of 
the incident radiation lie in the range commonly called the 
resonance region. The classical method of solution in the 
resonance region utilizes the separation of variables tech- 
nique. Following this approach Asano and Yamamoto [2] have 
solved the problem for spheroidal particles. Generalization 
of their results for the case of randomly oriented spheroidal 
particles is presented in [3,4]. 

Based on the well-known T-matrix approach, an integral 
equation method which was introduced by Waterman [4,5], 
Barber and Yeh [7] have solved the electromagnetic scattering 
problem by arbitrarily shaped dielectric bodies. A review 
article for scattering by non-spherical particles is listed 
for reference [12] . 

The scattering of microwaves by raindrops is also a 
scattering problem which was investigated by many re- 
searchers. It is clearly related to the problem of investi- 
gating and modeling the shape of raindrops (see Pruppacher- 
Klett [34] for extensive discussion) . The results for spher- 
ical raindrops and hail particles are given in Battan [8] . 
For the results which use the spheroidal description of rain- 
drops refer to Doviak and Zrnic [46] . Also, Warner and Hizal 
[44] have investigated the scattering of microwaves by 
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spheroidal raindrops based on an integral equation method. 
Recently Beard and Als [10,13] have proposed a simple model 
for the electrostatic equilibrium shape of falling raindrops. 
We also mention Jameson’s work [21,22,23] for radar measure- 
ments of rainfall and for estimation of raindrop size distri- 
bution. A review of the reflectivity technique of measuring 

rainfall is presented in [42] . 

In this work the scattering by ellipsoidal raindrops in 
low-frequencies is examined. At frequencies below 6 GHz (or 
wavelength \=5 cm) most of the rain droplet sizes satisfy the 
condition ka«l (k is the wave number and a is the character- 
istic diameter of the scatterer) and therefore Rayleigh scat- 
tering is applicable. The approximate upper limit of the 
characteristic radius of the scatterer is generally taken to 
be a = 0.05X [20]. At this radius the error of Rayleigh ap- 
proximation is less than 4% [27] . We assume that the rain- 

drop is an ellipsoid, with semi -axes a lf a 2 , and a 3 . Thus, we 
have to solve an electromagnetic scattering problem in SR 3 . 
Compared to the already existing solutions, we have one more 
degree of freedom. We assume an arbitrary direction of the 
incident radiation on the scatterer, and we examine the far- 
field patterns, taking into account the particle size distri- 
bution of the rain drops and their random orientation. All 
our analytical results are graphically compared to the known 
results for corresponding spheres and spheroids . 
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II. 2. 


MATHEMATICAL FORMULATION OF THE PROBLEM 


II. 2.1. Equations of electrodynamics and the fundamental 
dyadic solution 


We consider the propagation of electromagnetic waves in 
a medium. As it is well known, the electric field E(r,t) and 
the magnetic field H(r,t) are governed in free charge and 
current space by the Maxwell's equations 


VxE(r,t)-|i 


8H (r,t) 


3t 


(la) 


V *E (r, t) =0 


(lb) 


8E(r,t) _ 

V *H (r, t) =-e + oE(r,t) 

dt 


(2a) 


V-H(r,t)=0 

where e is the dielectric constant, p. the permeability and O 

the conductivity of the medium. 

For a general consideration of the electromagnetic prob- 
lem we refer to [34,39]. Assuming, without any loss of gen- 
erality, harmonic time dependence for the electric and the 
magnetic field, we have 

E(r,t)=E(r)e“ icot (3) 
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H (r,t)=H(r)e" itot 


(4) 


where CO is the angular frequency. In what follows, we can 
suppress the time dependence from all field quantities and 
then, for steady-state waves the equations corresponding to 
equations (1) and (2) are 


V xE (r) =i|lH (r) 


(5a) 


V -E(r)=0 


(5b) 


V xH(r) = (-eiCO + a) E (r) 


(6a) 


V *H (r) =0 


(6b) 


Elimination of the H(r) field in the equations (5) by substi- 
tution of the equations (6) gives us the following equations 
for the electric field 

V x V xE (r) -k^E (r) =0 (7a) 

V *E (r) =0 (7b) 


where k is the complex propagation constant 
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( 8 ) 


Jc2=o) 2 LIE(1 + i ) 
ecu 


In particular for nonconducting media (a - 0) k is real and 
it is expressed in terms of the phase velocity c as 

k = 0)(|i£)" 1/2 = 0)/c 

Similarly, the elimination of the E(r) field in the equations 
(6) gives us the equations of the same type as equations (7) 
for the magnetic field. 

The fundamental dyadic solut ion T (r, r ' ) satisfies the 
equation 

VxVxf (r,r')-k 2 f (r,r')— 47tl5(r-r' ) (10) 


where r is the position of the observation point, r’is the 
position of the source point, I is the identity dyadic, and 
S(r-r') is the three-dimensional delta function. After some 

derivations we conclude that 


f <r,r') 


e ik | r-r' I 
k 2 | r-r’ | 3 


{ k 2 (r-r' ) ® (r-r ' ) 


_ (r-r' ) <8> (r-r* ) 

+ (1— ik | r— r' I) [I 3 | r— it ' I 2 


_ e ik | r-r ' I .. 
_I I r-r 'I * 


(ID 
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IX. 2. FORMULATION OF THE SCATTERING PROBLEM 


Let us assume that V 2 is a bounded convex and closed 
subset of having a smooth boundary S. Let V 2 be a dielec- 
tric with dielectric constant £2 and permeability H 2 that lies 
in an infinite, homogeneous isotropic medium Vj_ with dielec- 
tric constant £i and permeability |li . We assume, as previ- 
ously, harmonic time dependence. An incident plane electric 
wave E in propagates in the medium V 1 along the propagation 

vector k. Let the corresponding magnetic wave be H in . The 
two waves have the form 

gin( r ) = jj e ikik*r (12) 

. ... A 

H in (r) * kxbVei/M-i e ik i k * r (13) 

where b is the unit polarization vector for the electric 
field so that b-k = 0 and k x is the propagation constant for 

Vi- 

If E(r), H(r) are the scattered electric and magnetic 
waves, respectively, and Ei(r), Hi(r) the total fields for 
the spaces Vi, i = 1,2, then due to linearity the total waves 
are given by the sum of the incident plus the scattered 
field. 

The vector fields Ei n (r), E (r) , Ej. (r) , Hi n (r), H(r), 
Hi (r) , satisfy the equations 
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(14) 


V x V x w (r) -k i w (r) =0 


reVii=l, 2 


V *w (r) =0 

where 

ki=O) 2 0iM.i (16) 

The boundary conditions for the electric field on the 
surface of the dielectric are given by the equation 

n xBi(r') = n xB 2 (r') r’eS (16a) 

n x [VxKi(r') ] = x[Vxii(r')] r'eS (16b) 

M-2 

On the surface of the scatterer the boundary condition 

n ‘Ei(r') = n *E2(r') r’eS (16c) 


must also be satisfied as a consequence of the integral rela 


tion (see [37]) 


| n • E (r ' ) dS (r ' ) = 0 
S 


(16d) 


49 


The scattered fields E(r), H(r) satisfy the radiation 


condition, due to Sommerfeld [27] : 

lim r x ( Vx ]) + ik i r [H(r)] = 0 (17) 

r— 

uniformly over all directions. 

For a general consideration of the electromagnetic scat- 
tering problem we also refer to [27] and [43] . 

The total electric field admits the following integral 
representation [1] 

Ei(r') - Ei*Mr) - — f [ (“ - 1) k* E 2 (r • ) • f (r, r ' ) 

4tc J ei 

v 2 


+ (l - — ) [V xE 2 (r ' ) ] V r • x T(r, r' ) ] dU (r ' ) (18) 

H2 

where the index jr 1 means differentiation with respect to the 
variable r. It also has been proved [28] that the normalized 
spherical scattering amplitude is given by the relation 


g(r,k) =— (iki)3[ J <— - 1) E 2 (r ' ) e - ik l r * r ' dU (r • ) • (I-r®* ) 

47C ' €l 

v 2 
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+ k. 


J <1 

v 2 


— )VxK2(r')e ikir*r' 
M-2 


dU(r' ) • I xr 


(19) 


where the normalized scattering amplitude is defined by the 
relation 

S ( r) = g (r, k) h(kir) + 0 (20) 


where h(kir) is the zeroth order spherical Hankel function of 
the first kind 


e ikir 

h(k l r) " TklT 


( 21 ) 


In radar applications, the bistatic radar cross section 
CT b i and the back-scattering cross-section a b are often used. 
They are related to normalized scattering amplitude through 
the relations 
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0bi = k2 


A A O 

lg(r,k) I 


( 22 ) 


CJ b 


47t 

k2 


/S A O 

lg(-k,k) r 


(23) 


The back-scattering cross section O b is also called the radar 
cross section. We define as scattering cross-section the ra- 
tio of the time average rate (over a period) at which energy 
is scattered by the body, to the corresponding time average 
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rate at which the energy of the incident wave crosses a unit 
area normal to the direction of propagation. The scattering 
cross-section is related to the normalized scattering ampli- 
tude via the relation 

0b .i J i,<-k,k)i 2 atm < 24 > 

1 * 1-1 


Xi. 3. THE PROBLEM IN LOW-FREQUENCIES 

It is possibile to work in low-frequencies when ka-0 
(i.e., 27caA=0) , where k is the wavenumber, X is the wave- 

length and a is the "characteristic dimension" of the scat- 
tered that is the radius of the smallest sphere that con- 
tains the scatterer. In such a case we can use the potential 
theory and approximate the problem by a sequence of potential 
problems. The potential theory approximation is also called 
long wavelength approximation or low-frequency approximation. 
The solutions of the vector Helmholtz equation considered as 
functions of the wavenumber are analytic in the neighborhood 
of zero. Thus, we can expand them in a convergent power se- 
ries, which we call low-frequency expansion. 

For the electric fields Ej.(r) for i-1,2 we have the ex- 
pansions 
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(25) 


■i(r) 


OO 


I 


(ikj) n 
n ! 



(r) 


Inserting this expansion into equation (14) and equating 
equal powers of k, the following sequence of partial differ- 
ential equations is obtained 

VxVxiji^ (r) + n(n-l)mi <)>n-2 = ® 


V * <t> < n > < r > = 0 


for n=0, 1,2,... 


(26) 


where 


mi = 


{m2 

Jiiei 


for i=l 
for i=2 


(27) 


The boundary conditions could be transformed into the bound 


ary conditions 


^lx(|> < J ) (r , ) - nx^ 1 ^’ (r' ) r'eS (28a) 

Sx[Vx* (1> (f)l - ^ Sx[Vx* < i > (r')l r' e S (28b) 

T n i » 0 11 
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r'eS 


(28c) 


The incident wave can also be expanded into a convergent 
power series of ki as follows: 


K in (r) = (i nt ) (*' r)n 

n-0 


The fundamental dyadic solution r(r,r') has the expansion 


- £ (iki) n ~ , 

r (r, r ' ) = b£ — ^7 — Yn< r ' r > 
n-0 


where 


y n (r,r’) = - 


' i n_1 r ~ (r-r' )® (r-r' ) t 

^ [ (n+1) I- (n-1) ——2 J (3D 


and the dyadic V r «xr(r, r') has the expansion 


~ v (iki) n £ , 

V r ' x r (r, r 1 ) = X n t ~ S n (r,r') 

n-0 


where 


S n (r,r') = (n-l)|r-r'| n ^(r,r') 
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Substituting (25), (29), (30), and (32) into (18) and equat- 

ing equal powers of ki the following integral relations among 
the coefficients (r) , . . . ,<|>^ n (r) are obtained 


( 1 ) 
1 n 


(r) = b* (k*r) n 


ft J [ (? - 0 p'p- 1 * *p-2 W 1 ' 1 ' 


p-0 v 2 £l 


^rr.^ ( 2) (rI) . Y n -p(r,r’) ] dU(r') (34) 
H2 p 


- (1 - — )Vx(t,' ' o 


In order to derive the low frequency expansion for the 

A A f 

scattering amplitude g(r, k) we need the expansion 

e ikir*r ' * £ (-i)ft (r-r') n < 35 ) 

n-0 


Substituting into equation (19) we obtain 


A A 

g(r,k) 


_i_ v < ili i )n ' 1 ' 3 ( n ) - i) 

« n t 0 n! 'P 7 


r (2) a ~ A A 

I 6 (-l)P (r.r')P dU(r' ) • (I-r®r) 

J T n-p 

v 2 
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_ 1 _ 

47C 


(p 

M-2 n =o P“° 


* (?) A — A 

f Vx* (r') (-l)P (r.r’)P dU(t')-lx* ( 36 ) 

J T n-p 

V 2 


In particular the leading term approximation as ki->0 is 

g(c,k) - - X (iki) 3 { J (- - l) 4>‘o > <r')dO(r')-(J-r®i) 
« V 2 61 

- f (l - — ) Vxp 1 ^ (r'ldO(r’)-iir } 

.{ V-2 

v 2 

+ 0<k*> < 37 > 

The leading term approximation for the scattering cross-sec- 
tion is given by 


a - - “ { (~ - l : 

6k ei 

) 2 1 j (r')dU(r’) 

V 2 

CM 

V-2 

J Vx<|) ( Q ) (r')dU(r') 1 2 } + 

v 2 

O(kJ) (38) 


In low— frequency regions the magnetic *field assumes the se- 
ries expansions: 
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re Vi 


1 = 1 / 2 


(39) 


H ± (r) 


oo 


s 


(iki ) n (i) 
n! V n 


(r) 


Following the same procedure as for the electric field 
we can arrive at a similar sequence of potential problems . 
We also mention the relation between the coefficients of the 
low-frequency series expansions for the electric and the mag- 
netic fields given by equations (25) and (39) , respectively 

Vx<b (l) (r) = (mi — ) 1/2 n \lf^!(r) reVi i=l,2 (40) 
n 


In Rayleigh scattering the radiation zone is determined 
from the relation 


S 


d£ 

X 


(41) 


where d is the characteristic diameter and X is the wave- 
length. So, we can assume a single scattering process if ev- 
ery scatterer lies in the radiation zone of any other scat- 
terer, that is at low-frequencies a distance equal to a cou- 
ple of diameters away from the scatterer justifies the appli- 
cability of our method. 
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II. 4. LOW-FREQUENCY SCATTERING BY AN ELLIPSOIDAL DIELECTRIC 

Let us assume that the triaxial ellipsoid 
2 

3 x. 

y — — <, l with 0<a3<a2 < ai <o ° (42) 

i-i a. 

is the dielectric scatterer. In order to reflect the geomet- 
rical peculiarities of the scatterer we introduce the ellip- 
soidal harmonic functions. 

II. 4.1. Ellipsoidal harmonic functions 


The ellipsoidal harmonic functions as it is well known, 
form a complete system of eigenfunctions. In what follows, 
we will give certain definitions about ellipsoidal harmonics. 
For details about the ellipsoidal harmonics we refer the 
reader to Hobson [19]. For details about the solution of the 
Laplace equation we refer to Morse and Feshbach [32] and for 
tabulated information of all the coordinate systems to the 

Moore and Spencer handbook [31] . 

The ellipsoidal coordinates (p,|l,v) are related to the 

Cartesian coordinates (x 1 ,X 2 ,X 3 > by 

v = WL (43a) 

1 ^2h3 
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(43b) 


P 2 - h 3 Aj 

f Jl 2 — h 2 

^ 3 -V 2 

x 2 = 

] 

fp 2 -h 2 

*l h 3 

^2' v2 


x 3 = h x h 2 


where 


h 


2 

1 



2 

a 3 


h 


2 

2 



2 

a 3 


2 2 2 
h 3 - a i - a 2 


(43c) 


(44a) 


(44b) 


(44c) 


and 


0^v 2 ^h3^)i 2 ^ h 2^P 2<00 


(44d) 


Separation of variables for the Laplace equation in ellip 
soidal coordinates produces the interior ellipsoidal harmon 
ics 


En(P.P»V) 


e”(P> e;(»i) E*"(V) 


m 


m 


( 45 ) 


and the exterior ellipsoidal harmonics 
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f”( P ,h,v) = f>) e“(ji) e>) 


(46) 


where E™ are the Lame functions of the first kind and 


F>) = ( 2n+l) E>) i”<P) 


(47) 


with 


oo 


m f 

du 


X n<P> = J 

/ 2 

/ 2 

p [E^(u) ] 2 ^ 

y u2 - h 2 ^ 

\j u2 ~ h 3 


(48) 


are the Lam6 functions of the second kind. The index n spec 
ifies the degree of the corresponding ellipsoidal harmonic 
and takes the value of n = 0,1, 2, 3,... while m represents the 
number of independent harmonic functions of degree n and runs 
through the values m - l, 2, . . . , 2n+l . In the present work we 
use the interior ellipsoidal harmonics of degree 0,1 and for 
the sake of completeness we give their exact form, both in 
ellipsoidal as well as in Cartesian representation: 


eJ (p,|l,V) = ppv = x L h 2 h 3 


(50a) 
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E^(p,M.,v) = ‘yj p 2 -^ ‘yj M- 2-h 3 *\J h 3 - ^ 2 “ x 2 h l h 3 (50b) 

E^(p,p,V) = p 2 "h2 yj^~\L 2 *\j 4^ 2 = x 3 h l h 2 (50c) 

The exterior ellipsoidal harmonics of degree 0,1 are given 

from equation (46) when equations (49) -(50) are used. The 

Lam6 functions of degree 0,1 that appear in the expression 
(48) for the elliptic integrals I n (p) are 

Eq (P> - 1 

E?(P> = 

The set 

{E >) E>): n=0, 1,2,... m=*l, 2, . . . , 2n+l } 

forms a complete orthogonal set of surface harmonics on the 
surface of the ellipsoid. 

II. 4. 2. The elliptic integrals 

The four elliptic integrals lQ(a x ), I^(a x ) for m=l,2,3 
given by equation (48) for p— a^ are related by the formula 


Y 


2 

V 


> 2 -a , +a, 


for m=l,2,3 


(51) 
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I 


i-1 


m 


Ii<ai) 


1 

a l a 2 a 3 


(52) 


t, 4 - ij<»i> 


(53) 


For other useful formulas related to elliptic integrals we 
refer to [17] . In order to express the elliptic integral 

lj(a!) to its canonical form we apply some transformations and 

conclude in notation of elliptic integrals of the first kind 
that 


Iq ( a i) 


sin<|>o 



dt 

^ l-t 2 sin 2 a 0 


(54) 


where 


Oq = sin -1 



2 2 
a l -a 3 

(55a) 

2 

a l 
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(55b) 



For 1^ (a^) , for example, we take in standard notation of el 
liptic integrals of the second kind the relations 


Iq (*i) 




-sin 2 ao 


[E (<)>o/ ao) -F Wo* a 0> 1 


(56) 


a l a 3 


2 3 

Front equations (52—56) we can evaluate I^(a^) and I^(a^) 


II . 4 . 3 • The zeroth order approximation for the electric 
field 

The zeroth order approximation for the electric field is 
the solution of the boundary value problem 

VxVx<(> ( Q ) (r) - 0 for i-1,2 (57a) 

V*<|>*Q*(r) = 0 for i=l, 2 (57b) 

nx<|> ( J > <r') - nx^ ( J ) <r») r’eS (57c) 
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r'eS 


(57d) 


nx[Vx^P ( r ' ) ] - — nx[Vx ( [) ( Q ) (r‘) ] 
u M-2 

n-^Q* (r') = ^- n *()) ( q ) (r’ ) r’eS 



(r) 


+ O(-) 
r 


(57f) 


If we use the well-known representation for the electrostatic 
problem and that a particular solution for the exterior field 

A 

is equal to b, we have 


♦ ^ (r) = b + VU ( J’ (r) 

<t> ( o > < r > * v ° ( o ) (r) 


(58) 


(59) 


where the scalar potentials for the exterior and the interior 
fields are given in terms of second and first kind ellip- 
soidal harmonics as follows: 


U 


(1) 

0 


(r) 


a 00 VP> 


+ Z a 0i* Fi (P*M-*V) 
m—1 


(60a) 


U 


( 2 ) 

0 


(r) 




(60b) 
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For the evaluation of V F™ (p,|l,V) we will use the general form 


VF®(p,|l,V) - (2n+l)VE”(p,^V)lJ(p) 


- (2n+l) 


_e_ 

hp 


E*(p,P,V) 


[E„(P) ] 2 *\j P 2 " h 2 'Sj P 2_h 


(61) 


So, 


* (1 o <*> 


(l)m 

v apl 

b + 3h 1 h 2 h 3 X hm I l ( P )x m 
m-1 


(l)m 
3 3-qi 

[C 1 * 3E — E>,E>>] <62) 


> i u a 00 '“ F m 

■y p2_(j,2 -yp2_v2 m-1 Ej^ (p) 


♦ ‘S’ <r> 


- *1*2*3 X 


. (2)m 

3 b 01 


m-1 


hm Xm 


(63) 


Applying the boundary conditions on p = ai (the surface 
of the ellipsoid) and by the orthogonality of the surface el- 
lipsoidal harmonics we obtain a system from which we can 

. (1)1 , (Dm j =1 o 

evaluate the unknown coefficients a 00 and a 01 tor i i, 

So, the zeroth order coefficients for the electric field are: 
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♦ ‘J’w -b 


-I [ 


km^m 


a l a 2 a 3 


( a -0 


£i 


m* 


' L 3h 1 h 2 h 3 /e 2 \ m 

1 a l a 2 a 3 \ “ x ) Ii( a l) +1 


V 


fJ (p,p,V) 


] 


( 64 ) 


♦ 


( 2 ) 

0 


(r) 



m*l 


hih 2 h3 

a l a 2 a 3 


1 

- l ) 

El 


V 


E® (P41.V) 


] 


(65) 


The zeroth order approximation for the magnetic 
field 


The zeroth order approximation for the magnetic field is 
the solution of a boundary value problem similar to that de- 
scribed by equation (57) . This is due to the invariance of 
the boundary conditions for the dielectric under the substi- 
tution So, the zeroth order coefficients for the mag- 
netic field are: 


V ( J ) (r) = kxb 


t ^ r (kxb) -xh^ei/pi 

"VEi/pi - 2 , L hih 2 h 3 


m-l 


66 


( 66 ) 


(U _ \ 

a l a 2 a 3 V 1 / 

M’l 


a l a 2 a 3 


(— - i) iT( a i)+i 
hi 


V F™ (p,|l,V) ] 


¥ < 0 > < r > 


Ji, r (kxb) >xh m Vei/^i 

L L h 1 h 2 h 3 


m-1 


a l a 2 a 3 “ 0 Ii( a i) +1 


V E™(p4l,V) ] 


(67) 


The leading term approximation for the normalized 
scattering amplitude 


In order to evaluate the leading term approximation for 
the normalized scattering amplitude, given by Equation (37) 
we have to evaluate the integrals which appear in that equa- 
tion. First, we exploit the relation between the electric 
and the magnetic low-frequency coefficients given by Equation 
(40) . We have 


Vx<t> ( ^ (r) 



( 2 ) 

0 


(r) 


( 68 ) 


From the zeroth order approximation for the electric and 
magnetic fields we obtain 
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( 69 ) 


f ( 2 ) 471 

j V q J (r’)dU(r') = 3- a!a 2 a3 

V 2 


X b m*m 
m«l 


1 

a l a 2 a 3 (“ ~ l) I l< a l) + 1 

El 


J Vx<|> ( ^ ) (r')dU(r’) 

v 2 


4 31 JX2 

— a]a 2 a 3 




(kxb) -Xm 


i (70) 

a l a 2 a 3 (“ " 1 ) I l (a l )+1 


The normalized scattering amplitude is given by the re 
lation 


A A 

g(r,k) 


- (iki) 3 


a l a 2 a 3 

3 





a l a 2 a 3 


(52 - l) 

8 l 


(I-r®r) 
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1 


+ 




(Jcxb) -x ra 


a l a 2 a 3 (jj^ “ 0 Ii( a l) +1 


A A 

Ixr 


] 


+ O (k, ) 


( 71 ) 


II. 4. 6. The leading term approximation for the scattering 


cross-section 


For the scattering cross-section the following can be 


obtained 


_ 871 ( 1 ^) 4 
CT s “ 27 



2 2 2 
a l a 2 a 3 


[a l32 a 3 (^- l)l>i>*l ] 2 

El 

+ (l-^) 2 a l a 2 a 3 1 

111 m-1 

i }+0(k*) (72) 

[ a l a 2 a 3 ({“ " 0 < a l> +1 ] 

ri 
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II. 4. 7. 


The back-scattering cross-section 


Substituting in Equation ( 23 ) we conclude that the lead 
ing term approximation for the radar cross-section is given 
by the relation 


8Jt(k!) 4 2 2 2 V 

Ob = 9 a l a 2 a 3 2* 




- k_ A-k 1 + Bxk 


v m 


0<k®) 


( 73 ) 


where 


A = (A lf A 2 f A 3 ) and B = (B 1 ,B 2 ,B 3 ) 


(74) 




<?-i) 

El 




a l a 2 a 3 


(- - i) iT< a i)+i 

El 


( 75 ) 


B m = 


/M£ _ 

' 111 


17 (kxb) ‘Xjjj 


a l a 2 a 3 ” l) Ix(ai) + 1 


( 76 ) 


Two cases which are of special interest are examined in 
the sequel. The first is the radar scattering cross-section 
with vertical polarization, that means 
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and the second case of the horizontal polarization is 


(kxb) **3 = 0 


( 78 ) 


In these cases the leading term of the backscattering cross 
section is also given by Equation (73), but for the vertical 
polarization we have 


and for the horizontal polarization we conclude from Equation 
(78) that 


B 3 = 0 


(80) 


II. 5. PARTICLE SIZE DISTRIBUTION 

Up to this point we have examined the scattering by a 
single ellipsoid. However, our main interest lies in how a 
wave interacts with many randomly distributed particles. The 
particles we deal with in practice are not usually all one 
size but normally their sizes are distributed over a certain 
range. It is, therefore, important to take into account the 

size distribution of the particles. 

Let n(D)dD be the number of particles per unit volume 
having a dimension (such as diameter) between D and D+dD. 
The total number of particles per unit volume is then 
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( 81 ) 


P “ 


OO 



which is called number density (or simply density) . 

We can also define by w(D) a probability density func- 
tion for finding the particle size between D and D+dD 


w (D) 


n (D) 

P 


( 82 ) 


where 


OO 

J w(d)dD = 1 
0 


( 83 ) 


Now we can define the average cross-section and the av 
erage radar cross-section as the following 


OO 

E{a} = J C(D)w(d)dD 
0 


( 84 ) 


and 


E{a b } 


OO 

J O b (D)w(d)dD 
0 


( 85 ) 
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where E{} is the expectation operator and a, a b are given by 

Equations (72) and (73), respectively. These quantities can 
be expressed in terms of the semi-axis of the ellipsoid a lf 

ci and a^ • 

The size distribution n(D) in Equation (81) can be rep- 
resented by an exponential distribution or by a three parame- 
ter gamma distribution. 

In order to take into account the particle size distri- 
bution, using equations (84) and (85), we need to establish a 
relation between the semi -axes of the ellipsoid (in terms of 
which are expressed the scattering cross-section and the 
radar cross-section) and the radius of an equivalent volume 
spherical raindrop, because in terms of this parameter we 
have the information on the particle size distribution. 

The principal curvatures of a point on the surface of an 
ellipsoid p=a x are given by the relations 


a l a 2 a 3 


a i - M- 2 y a i~ v2 


*-ii2 


( 86 ) 


= a l a 2 a 3 1 — (87) 

2 " a? - v2 

The sum of the principal curvatures at the points on the 
equator of the ellipsoid are given by the form 
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( 88 ) 


+ k 2 - a 2 a 2 [a 2 (a 2 + a 2 ) + <a 2 - ajlxj 
(a: 2 [a? ♦ (a 2 - a 2 >x 2 ]- 2 ' 2 } 


For an oblate spheroid a.\ - a 2 and 

a i 1 

ki + k 2 = + “ 

a 3 a 2 

which is the same result as in the Appendix of [ 18 ] . 

= 0 i.e. at the point (0,a2,0) we have 

ki + k 2 ** a2 (“^2 + “ 2 ) 
a 3 a l 

and for x x - a x (the same as for x 1 - - ai), that is 
point (a lf 0,0) f we obtain 

kl + k 2 = ai (“^2 + ~ 2 ^ 
a 3 a 2 

It holds that 

a 2 S + k 2 £ ai (-7 + -7) 

a 3 a x a 3 a 2 

If we assume that 


( 89 ) 


For xi 


( 90 ) 


at the 


( 91 ) 


( 91 ) 
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(92) 


x 

a i 


where parameter Xe(0,l), we can study the influence of defor 
mation of the spheroid along the 

We will follow the simplified analysis of the raindrop 
shape problem due to Green [18] . From the mechanical equi- 
librium condition on the surface of the drop we have 


where R L , R 2 are the principal radii of curvature, and Pi and 
P e are internal and external pressures, respectively. 

By ignoring the aerodynamic and hydrodynamic pressures 
at the equator of the spheroid Green [18] established the 
following relation 

a (— + “) = 2aaQ 1 + pg'b = a(ab“ 2 + a -1 ) (94) 


where a Q is the radius of the sphere with equal volume as the 
oblate spheroid with axes a 1 =a 2 =‘a, and b=a 3 . The right hand 
side of Equation (94) is constant and independent of the 
point on the equator of the spheroid due to the fact that the 
equator is a circle. At the points of the equator of the el- 
lipsoid the sum 
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( 95 ) 


R 1 


+ 


_ 1 _ 

r 2 


k L + k 2 


which is obvious from Equation ( 91 ) and depends on the spe- 
cific point. Therefore, in order to establish a relation be- 
tween the semi-axes a 1 /a 2 ,a 3 and the radius of the equivolume 

sphere based on the physics of the problem, we first choose 
the "mean— value" of the sum of curvatures 


||>1 + k 2 > I (ax ,o, 0 > + (k l + k2) ' ( 0 , a 2 , 0 , 


2 2 2. -l 


-o') 23 

[ a^a^ (a x + a 2 ) + a 3 (a 1 + a 2 ) J (2a 1 a 2 a 3 ) 


( 96 ) 


Now we can establish the equality of the volumes of the el 
lipsoid and the sphere, through the relation 

3 (97) 

a l a 2 a 3 = a 0 


So, with the same arguments as Green we conclude that 



2 3 

+ 33(3! 


+ a 2 ) ] ( 23 ^ 2 ^ 3 ) 


2+Br 1 ( 98 ) 

a 0 


where 


B = pg’ 


2 

a 0 


a -1 


( 99 ) 
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is the Bond number . 

If we introduce a new variable 

5 - (2i) 2 - 1 < 100 > 

a 0 

we have that 

ai = aoVTTTiT < 101) 

a 2 = Xa 0 V (1+5) (102) 

a 3 = X-iao V (1 + S)~ (103) 

Substituting in Equation (98) Equations (8), (15), (18), 

(96), (100), and (103) we obtain 

B - Vu7«r[ xi< v X) ' « i + 5) + i ir ] - 2X(1 + 5) 

(104) 

If we take that 

■\j (l + 8) = 1 + j + 0(8 3 ) (105) 

we conclude that Equation (104) can be represented to 0(8 2 ) by 
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.. 35X 4 (1 + X) - (1 + X 3 ) 

B ’ 53 To 

7X 4 (1 + X) + (1 + X 3 ) - 8X 2 

+ 5 ~a 

X 4 (l + X) + (l + x 3 ) - 4X 2 (106) 

+ 2X 

For X=1 (the case of the spheroid) the same relation as 
Green's is obtained. From Equation (104) we can evaluate 5 

in terms of the Bond number and from Equations (101), (102), 

and (103) we have a 1 ,a 2 ,a 3 in terms of the radius of the 

sphere with equal volume. Thus, from Equations (84) and (85) 
we can evaluate the average cross-section and the average 
backscatt ering cross-section, respectively. 


II. 6 . RANDOMLY ORIENTED ELLIPSOIDAL PARTICLES 

In order to take into account the orientation of the 
scattering particles we will take the average depending on 
the orientation. Thus, we choose a reference rectangular 
Cartesian coordinate system and introduce as un_ 
known the Euler angles of the transformation (by rotation) 
from the reference system to one coinciding with the princi- 
pal axes of the ellipsoid. For the unit vectors of the 
(Yi f y 2 , y 3 ) system in terms of x L (the unit vectors of the sys- 
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tem which coincides with the principal axes of the ellipsoid) 


we have 


y - 


X'D 


for i = 1/2,3 


(107) 


where the elements of the matrix D are functions of the Euler 


angles and are given by the relations 


d 13 ^ = sin©! sin0 2 + 

COS0! 

COS0 2 

COS0 3 

di2 * cos©! sin0 2 ~ 

sin0! 

COS0 2 

COS0 3 

d i3 = cos0 2 sin0 3 




d 2 i ” sin0! cos0 2 + 

COS0! 

sin0 2 

COS0 3 

d 22 = COS0! COS0 2 + 

sin0 x 

sin0 2 

cos© 3 

d 2 3 =* " sin0 2 sin0 3 





d 31 = - cos©! sin0 3 


d 32 * sin0i sin0 3 


d 33 = cos0 3 


(108a) 

(108b) 

(108c) 

( 108d) 
(108e) 

( 108f ) 
( 108g) 
( 108h) 
( 108i) 


where O^Q^n/2 for i=l,2,3. 
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We mentioned before that the angle 0 3 is the angle be- 
tween the y 3 and the x 3 axes. The x x x 2 plane intersects the 
y lY 2 plane in a line, which is called the nodal line. The 
angle is the angle between the nodal line and the y^-axis and 
q 2 is the angle between the x^axis and the nodal line. 

The relations between the components of the vectors rel- 
ative to the reference frame and the frame coinciding with 
the principal axis of the ellipsoid are 

k* = X cljmkj ^ m = 1,2,3 (109) 

m— 1 


km 


^ 1 

X j 


for m = 1,2,3 


( 110 ) 


where 


A 


k 


X ktn*m 


m-1 


X 

m-1 


1 A 1 

m 


( 111 ) 


b 


m-1 


3 A 3 • A » 

X = X k m y m 

m-1 


( 112 ) 


So, in order to take the scattering cross-section and the 
radar cross-section over all the orientations we obtain 
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<CT> 


( 113 ) 


JC/2 Jt/2 Jl/2 

= J j 1 O (0!, 0 2 , 63 ) f (9lf ®2r d0id02d03 

0 0 0 

where f( 0 lf 62 , 63 ) is the probability density function for the 
angles ( 61 , 62 , 63 ) and a(0 1 ,e 2 ,0 3 ) is the scattering cross-sec- 
tion given by Equation (72) after the substitution of the 
components k m , b m in terms of k m , b m . 

Similarly for the backscattering cross-section we have 
the average over orientation 

jt/2 jt/2 n/2 

«J b > = j | j O b (0 1 ,02,0 3 ) f (01 ^ 2 ^ 3 ) (114) 

0 0 0 ” 

If we want to obtain the average over the size and ori- 
entation we must substitute in Equations (113) and (114) the 
terms O and a b by the averages given by Equations (84) and 

(85), respectively. 

II. 7. NUMERICAL RESULTS - DISCUSSION 

The solutions presented in the above sections can be 
tested numerically against the known solutions for their spe- 
cial cases. For example, if ai=a 2 =a 3 we can use the solution 

of Rayleigh [ 8 ] . Figure 2 in the Appendix A presents both 
solution and the differences are indistinguishable. Another 
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spectial case of interest in radar-rainfall estimation is 
that of an oblate spheroid. Plots in Appendix A present the 
comparison with the Gans ’ solution for both single scatterer 
and volume scattering cases. Computer programs written in 
FORTRAN are included in Appendix B. 


82 



1 1. 8. FURTHER RESEARCH 


In the present work we have examined the scattering of 
raindrops in the Rayleigh region. We have assumed the shape 
of the large drops as an ellipsoidal one and we have calcu- 
lated the Rayleigh approximation up to the first order. In 
order to improve the accuracy of our results we have to take 
into account the second order approximation of the Rayleigh 
series as well. Further, for the particle size distribution, 
when we calculate the sum of the curvature, at any point of 
the ellipsoid, we have to assume a more ''realistic” average 

value of the sum of the curvatures. 

In order to derive numerical results taking into account 
the random orientation of the raindrops as well, we have to 
think about the probability density function, for the Euler 
angles. Beard and Jameson’s results about the canting of the 
raindrops give us enough information for the probability den- 
sity function for one of the three Euler angles [9] . 

If we want to have more accuracy for the scattering 
problem of the large raindrops, we have to assume that their 
shape is better approximated by a spherical cap. The only 
known results for spherical caps are due to Thomas [40] for 
the acoustical case and to Collins [14,15] for the electro- 
magnetic case. These last results have been obtained assum- 
ing different boundary conditions from that which we consider 
on the raindrop surface. From a mathematical point of view 


83 


such a study would be rather difficult, but the results would 
be applicable for the entire range of frequencies . 

The consideration of the scattering of the raindrops as 
a multi-scattering problem is also another possibility. Such 
an investigation can be based on Twersky's work [41]. We 
also can examine the possibility of exploiting the results of 
Peterson and Strom [33] who have generalized the T-matrix ap- 
proach of Waterman for the case of multiple scattering by N 
arbitrarily shaped configurations. 
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SCATTERING FROM WATER SPHERES 
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*** subroutine calculates the radar cross-section for an 
ellipsoidal scatterer (Kiriaki-Krau ewski method) 

(uses IMSL special functions) 


Parameters : 
am 
rk 
bm 

rmag_sus 
epsl 
eps2 
r lambda 
sigmab 


- vector of ellipsoid axes (m) 

- unit vector of incident wave direction 

- unit vector of polarization 

- magnetic suseptibility of the scatterer 

- dialectric constant for the scatterer 

- dielectric constant for the medium 

- wavelength (m) **-m 

- backscattering crossection (output, m 2 .) 


implicit real*8 (A-H, 0-Z) 
dimension A (3) ,am(3) ,B(3) , 


rk (3 ) , bm ( 3 ) ,RI1(3) ,h(3) 


pi=DCONST( ' PI | ) 
constant =4 . *pi/9 . 
wave_num= 2 . 0*pi/rlambda 


x=am(l) *am(l) 
y=am(2) *am(2) 
z=am(3 ) *am ( 3 ) 
al23=am(l) *am(2) *am(3) 

h(l) =sqrt (y-z) 
h (2) =sqrt (x-z) 
h (3 ) =sqrt (x-y) 

sss=DMACH(l) 

bbb=DMACH(2) 

RI1 (1) =DELRD (y , z , x) /3 .0 
RI1 (2) =DELRD(x, z,y) /3 .0 
RIl (3) =DELRD(x,y, z) /3 .0 

sum=RIl ( 1 ) +RI1 ( 2 ) +RI1 ( 3 ) 
zp=l . 0/ ( am ( 1 ) * am ( 2 ) *am(3) ) 

coef f l=eps27epsl -1. 
coef f 2=4 . 0*pi*rmag_sus 
coef f 3 =al23 *coef f 2 


do m=l, 3 

ar=bm(m) *coef fl 

A (m) =ar/ (al23 *coef f 1*RI1 (m) +1. ) 
enado 


★ 


★ 


★ 


B (1) =coef £ 2 * (rk (2 ) *bm(3 ) -rk(3) *bm(2) ) 
/ (coeff3*RH (1) +1. ) 

B(2) =coeff2* (rk(3) *bm(l) -rk(l) *bm(3) ) 
/ (coef f3*RIl (2) +1 . ) 

B ( 3 ) =coeff2*(rk(l) *bm(2) -rk(2) *bm(l) ) 
/ (coef f3*RIl (3 ) +1 . ) 


sum=0 . 0 
do i=l,3 

sum=sura+A ( i ) *rk ( i ) 
enddo 


sl= (A (1) -rk (1) *sa+rk (3 ) *B(2)-rk(2)*B(3))**2 
s2=(A(2)-rk(2) *sa+rk(l) *B(3) -rk(3) *B(1) ) **2 
s3= (A (3) -rk (3 ) *sa-rk(2) *B(1) -rk(l) *B(2) ) **2 


zp=x*y*z 

s igmab=wave_num* * 


4*constant*zp* (sl+s2+s3) 


return 

end 



n o o o 


c* 

c* 


★ 


★ 

c* 


c* 

c* 




★**★***★★**★**★** 


***★★★****************' 



COMPLEX crindx 


★ * * 


SUBROUTINE COMPUTES THE NORMALIZED CROSS SECTION FOR 
A SPHERE IN RALEIGH SCATTERING REGIME 


complex zp,K 
zp=crindx*crindx 

K= ( zp-1 . 0 ) / (zp+2.0) 
RK2=ABS (K) *ABS (K) 

sigma=4 . 0*alfa**4*RK2 

return 

end 





c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 


c 


:★★*★★★★**★*********** 


**★★*★★★★★★*********** 

subroutine bhmie (crindx, alfa, qext , ^sca^qbaciO 

***************************** 

complex crindx 


********************** 


************************************* 


Subroutine calculates Mie scattering by water spheres. 

It is based on the program included in Bohren & Hoffman (1983 


Parameters: . . , , . _ . . 

crindx ~ complex refractive index (relative) 
of the scatterer 

alfa - the size parameter (p*pi*radius/wavelength) 

qext - total extinction efficiency 

qsca - total scatering efficiency 

qback - normalized backscattering crossection 


dimension emu (100) , theta ( 100) . pai (100 l.tau (100 ) piO ( 100) , pil (100) 
complex d(3000) ,y,xi,xi0,xil, an, bn, si (200) ,s2(200) 
double precision psiO , psil , psi , dn, dx 


pi=CONST ( ' PI' ) 
dx=alf a 
y=alf a*crindx 

*** series terminated after nstop terms 


nang=ll 

xs top=alf a+4 . *alf a* *0.3333+2.0 

nstop=xstop 

ymod=cabs (y) 

nmx=amaxl (xstop,ymod) +15 
dang=pi/2 ./float (nang-1) 

do 100 j =1, nang 

theta ( j ) = ( float ( j ) -1 . ) *dang 
amu ( j ) =cos ( theta ( j ) ) 

100 continue 


*** logarithmic derivative d(j) calculated by downward 
recurrence beginning at j=nmx 

d ( nmx ) =cmp lx ( 0 . 0 , 0 . 0 ) 
nn=nmx-l 


c 

do 110 n=l , nn 
rn=nmx-n+l 

d(nmx-n) = (rn/y) - (1./ (d (nmx-n+1) +rn/y) ) 
110 continue 


c 


do 12 0 j =1 , nang 
piO ( j ) =0 . 0 
pil ( j ) =1 . 0 
120 continue 

nn=2*nang-l 

do 130 j=l, nn 

si ( j ) =cmplx ( 0 . 0 , 0 . 0 ) 
s2 ( j ) =cmplx (0 . 0 , 0 . 0 ) 

130 continue 

*** Riccati-Bessel functions with real argument x 
calculated by upward recurrence 

psi0=dcos (dx) 
psil=dsin (dx) 
chi0=-sin (alfa) 
chil=cos (alfa) 
apsi0=psi0 
apsil=psil 

xi0=cmplx (apsiO , -chiO ) 
xil=cmplx (apsil, -chil) 
qsca=0 . 0 
n=l 

135 continue 
dn=n 
rn=n 

fn= (2 . *rn+l . ) / (m* (rn+1. ) ) 
psi= (2 . *dn-l . ) *psil/dx-psi0 
apsi=psi 

chi= (2 . *rn-l. ) *chil/alfa-chi0 
xi=cmplx (apsi , -chi) 

an= (d (n) /crindx+rn/ alfa) *apsi-apsil 
an=an/ ( (d (n) / crindx+rn/ alfa) *xi-xil ) 
bn= (crindx*d (n) +rn/alfa) *apsi-apsil 
bn=bn/ ( (crindx*d(n) +rn/ alfa) *xi-xil) 

qsca=qsca+ (2 . *rn+l . ) * (cabs (an) *cabs (an) +cabs (bn) *cabs (bn) ) 


c 


140 


do 14 0 j =1 , nang 
jj=2*nang-j 

pai ( j ) =pil (j) . . 

tau ( j ) =rn*amu ( j ) *pai (j ) - ( rn+1 . ) piO ( d ) 


p= ( -1 . ) ** (n-1) . 

si ( j ) =sl ( j ) +fn* (an*pai ( j ) +bn*tau(j ) ) 


t = ( -1 . ) **n , . 

s2 ( j ) =s2 ( j ) + fn* (an*tau ( j ) +bn*pai (j ) ) 
if(j.eq.jj) go to 140 

si ( j j ) =sl ( j j ) +fn* (an*pai ( j ) *p+bn*tau( j ) *t) 
s2 ( j j ) *s2 ( j j ) +fn* (an*tau ( j ) *t+bn*pai ( j ) *p) 
continue 


psi0=psil 

psil=psi 

apsil=psil 

chi0=chil 

chil=chi 

xil=cmplx (apsil, -chil) 

n=n+l 

rn=n 


do 150 j=l, nang . 

pil ( j ) = ( (2 . *rn-l . ) / (rn-1 . ) ) *amu( j ) *pai ( 2 ) 
pil ( j ) =pil ( j ) -m*pi0 ( j ) / (rn-1. ) 
piO ( j ) =pai ( j ) 

150 continue 


if (n-l-nstop.lt .0) then 

go to 135 
else 

qsca= (2 . / (alfa*alfa) ) *qsca 

qext= (4 . / (alfa*alfa) ) *real (si (1) ) ... 

qback= (4 . / (alfa*alfa) ) *cabs (si (2*nang-l) ) *cabs (si (2 nang 1) ) 

endif 


return 

end 


APPENDIX C 


- BIAS simulation study 


92 




